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ABSTRACT 

The  radiative  transfer  (RT)  approach  is  widely  used  for  studying  scattering  from  layered  random  media  with 
rough  interfaces.  Although  it  has  been  successful  in  several  applications  such  as  remote  sensing  of  atmosphere 
it  is  well  known  that  this  approach  involves  certain  approximations.  In  this  paper  these  assumptions  and 
approximations  are  reexamined  and  explained.  To  enable  this  a  statistical  approach  is  employed  to  this  problem 
and  the  governing  equations  for  the  first  and  second  moments  of  the  wave  fields  are  derived.  A  transition  is 
hence  made  to  arrive  at  a  system  of  equations  corresponding  to  that  of  the  RT  approach.  It  is  thus  found  that 
more  conditions  are  implicitly  involved  in  the  RT  approach  than  generally  believed  to  be  necessary. 

Keywords:  random  media,  rough  surfaces,  radiative  transfer,  multilayer,  wave  approach 

1.  INTRODUCTION 

Atmosphere  is  often  modelled  as  a  stratified  medium  whose  scattering  and  absorption  characteristics  are  con¬ 
tinuous  function  of  height.  The  refractive  index  of  atmosphere  hence  has  a  continuously  varying  deterministic 
part  with  a  relatively  small  randomly  fluctuating  part.  In  addition  we  may  have  to  take  into  consideration 
reflection  from  the  irregular  ground.  In  this  connection  we  are  interested  in  studying  scattering  from  a  layered 
random  medium  with  rough  interfaces.  This  is  a  problem  that  one  encounters  often  in  many  applications  be¬ 
sides  atmospheric  sensing.  A  simple  approach  is  to  incoherently  add  the  contributions  of  volumetric  and  surface 
fluctuations. 1-3  However,  this  is  valid  only  when  we  are  in  the  single  scattering  regime.  There  are  some  other 
hybrid  approaches4  which  take  into  consideration  some  multiple  scattering  effects.  Brown5  outlines  an  iterative 
procedure  which  properly  includes  all  multiple  scattering  interactions.  However,  it  does  not  appear  feasible  to 
carry  out  the  calculation  beyond  one  or  two  iterations.  Among  the  other  methods  currently  used,  perhaps  the 
most  widely  used  approach  is  the  radiative  transfer  (RT)  approach. 6-9  Here  one  formulates  the  scattering  and 
propagation  in  each  layer  by  using  the  radiative  transfer  equation  which  involves  only  the  parameters  of  the 
medium  of  that  layer.  The  boundary  conditions  are  derived  separately  using  some  asymptotic  procedure  devel¬ 
oped  in  rough  surface  scattering  theory. 10-12  The  RT  equations,  along  with  the  boundary  conditions,  comprise 
the  system  that  describes  the  problem. 

Although  this  procedure  appears  to  be  reasonable  and  sound  it  is  apparent  that  certain  approximations  are 
involved  and  we  would  like  to  know  the  conditions  under  which  this  kind  of  approach  is  appropriate.  One  way 
to  better  understand  the  RT  approach  is  to  compare  it  with  the  more  rigorous  wave  approach.  For  the  case  of 
unbounded  random  media  it  was  found  that  the  RT  approach  is  applicable  when  (a)  one  uses  quasi-uniform  field 
approximation,  (b)  the  medium  is  statistically  quasi-homogeneous,  and  (c)  one  uses  the  ladder  approximation 
to  the  intensity  operator  of  the  Bethe-Salpeter  equation.13 

However,  our  problem  has  bounded  structures  which  are  randomly  rough.  Therefore  it  remains  to  be  seen 
whether  the  conditions  arrived  at  in  the  case  of  unbounded  random  media  will  be  sufficient  for  our  problem. 

In  this  paper  we  employ  a  wave  approach  using  surface  scattering  operators14  to  derive  the  transport  equations 
for  our  multilayer  problem.  In  this  process  we  find  that  there  are  more  conditions  implied  when  we  choose  to 
apply  the  RT  approach  to  our  problem  than  it  is  widely  believed  to  be  necessary.  One  such  condition  is  the  weak 
surface  correlation  approximation.  This  means  that  the  RT  approach  places  certain  restrictions  on  the  type  of 
rough  interfaces  that  it  can  model  accurately. 
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The  paper  is  organized  as  follows.  In  Section  II  we  describe  the  geometry  of  the  problem.  In  the  next  section 
we  describe  the  RT  approach  to  the  problem.  Next  we  describe  the  wave  approach  to  this  problem.  In  Section 
V  we  transition  to  the  transport  equation  system.  The  paper  concludes  with  a  discussion  of  our  findings. 

2.  DESCRIPTION  OF  THE  PROBLEM 

The  geometry  of  the  problem  is  shown  in  Figure  1.  We  have  an  iV-layer  random  medium  stack  with  rough 
interfaces.  The  permittivity  of  the  j-th  layer  is  tj  +  ej  (r)  where  €j  is  the  deterministic  part  and  ij  is  the  randomly 
fluctuating  part.  The  permeability  of  all  the  layers  is  that  of  free  space.  The  randomly  rough  interfaces  are  given 
as  z  =  Zj  +  (j( rj_).  It  is  assumed  that  and  (j  are  zero- mean  isotropic  stationary  random  processes  independent 
of  each  other.  Thus,  on  the  average  the  interfaces  are  parallel  planes.  Let  zq  =  0,  and  let  dj  be  the  thickness  of 
the  j-th  layer.  The  media  above  and  below  the  stack  are  homogeneous  with  parameters  eo,  ko,  and  ejv+i,  fciv+i, 
respectively.  This  system  is  excited  by  a  monochromatic  electromagnetic  plane  wave  and  we  are  interested  in 
formulating  the  resulting  multiple  scattering  process. 


f 

CN+ 1 


Figure  1.  Geometry  of  the  problem. 


3.  RADIATIVE  TRANSFER  APPROACH 

Multiple  scattering  in  a  complex  environment  is  well  described  by  the  radiative  transfer  theory.  This  theory  is 
not  only  conceptually  simple  but  also  very  efficient.  The  fundamental  quantity  here  is  the  specific  intensity  I 
which  is  governed  by  the  following  equation15-1' 

a  •  VI(r,  a)+  7  I(r,  a)  =  J  d£l'  P  (a,  a')I(r ,  s').  (1) 

One  may  regard  this  equation  as  a  statement  of  conservation  of  energy  density  I  which  is  a  phase-space  quantity  at 
position  r  and  direction  s.  7  is  the  extinction  matrix  which  is  a  measure  of  loss  of  energy  due  to  scattering  in  other 
directions.  P  is  the  phase  matrix  representing  increase  in  energy  density  due  to  scattering  from  neighbouring 
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elements,  fi  is  the  solid  angle  subtended  by  s.  Given  the  statistical  characteristics  of  the  medium  one  can  readily 
calculate  the  phase  function.  The  extinction  matrix  is  hence  calculated  using  the  optical  theorem.  The  specific 
intensity  in  each  layer  is  governed  by  an  equation  similar  to  (1).  Since  our  layer  problem  has  translational 
invariance  in  azimuth  the  RT  equation  for  the  m-th  layer  takes  the  following  form, 


COs9  —  Im(z,  s)+  7r 
dz 


I m(z,s)=  /  dfl'  Pm  (s,s')Im(z,s') 


(2) 


where  the  subscript  m  denotes  that  the  quantity  corresponds  to  those  of  the  m-th  layer  and  9  is  the  elevation 
angle  of  .§.  This  set  of  RT  equations  is  complemented  by  a  set  of  boundary  conditions  which  are  in  turn  based 
on  energy  conservation  considerations.  To  be  more  precise  we  impose  the  condition  that  the  energy  flux  density 
at  each  interface  is  conserved.  This  leads  to  the  following  boundary  conditions  on  the  m-th  interface 


I  2,(w) 


dfl  Rm+l,m  (^,  ^  ^  )  -p 


dfl  Tm,m+1  (^5  ^  $  )• 


(3) 


The  boundary  conditions  on  the  (m  —  l)-th  interface  are  given  as 


1)  s) 


dfl  R,m—  l,m  (£,  S  S  )  + 


dfl  Tm,m-1  (s,  5  1,  S  ), 


(4) 


where  Rmn  and  T mn  are  the  local  reflection  and  transmission  Mueller  matrices.  To  be  more  specific,  R mn 
represents  the  reflection  matrix  of  waves  incident  from  medium  n  on  the  interface  separating  medium  m  and 
medium  n.  The  superscripts  u  and  d  indicate  whether  the  intensity  corresponds  to  a  wave  travelling  upwards 
or  downwards.  The  integration  in  these  expressions  are  over  a  solid  angle  (hemisphere)  corresponding  to  s' . 
Suppose  we  have  a  plane  wave  incident  on  this  stack  from  above.  Then  the  downward  travelling  intensity  in 
Region  0  is 

Io(z,s)  =  B0<5  (cos 9q  -  cos 9Z)  d(</>0  -  </>i),  (5) 

where  Bo  is  the  intensity  of  the  incident  plane  wave  and  {9i,4>i]  describes  its  direction.  Since  there  is  no  source 
or  scatterer  in  Region  N  +  1, 

Iat+i(2>  s)  =  0. 

Notice  again  that  these  boundary  conditions  represent  conservation  of  intensity  at  the  interfaces.  We  should  point 
out  that  the  radiative  transfer  approach  as  applied  to  a  particular  problem  is  only  a  model  based  on  certain 
assumptions.  Since  the  RT  theory  is  used  in  a  variety  of  applications,  the  particular  assumptions  involved 
are  described  in  terms  of  different  terminologies,  specific  to  the  discipline  where  it  is  used.  One  good  way  to 
understand  in  more  general  terms  the  RT  approach  and  the  underlying  assumptions  is  to  compare  it  with  the 
more  rigorous  wave  approach.  For  the  case  of  an  unbounded  random  medium  this  kind  of  study  was  carried 
out  in  the  seventies.13  From  that  study  we  learn  that  the  radiative  transfer  theory  can  be  applied  under  the 
following  conditions: 

1.  Statistical  homogeneity  of  the  medium  fluctuations. 

2.  Quasi-stationary  field  approximation. 

3.  Weak  fluctuations: 

(a)  Ladder  approximation  to  the  intensity  operator. 

(b)  Kraichnan  approximation  to  the  mass  operator. 

These  are  the  well-known  conditions  that  we  associate  with  the  RT  approach.  However,  our  problem  has  bounded 
structures  and,  further,  they  are  randomly  rough.  The  question  is  this:  are  the  above  conditions  sufficient  to 
apply  the  RT  approach  for  our  problem?  This  is  the  motivation  for  this  paper.  We  follow  the  wave  approach  to 
this  problem,  derive  the  equations  for  the  intensities,  and  hence  make  the  transition  to  the  RT  equations.  This 
procedure  enables  us  to  better  understand  the  necessary  conditions  for  using  the  RT  approach  for  our  problem. 


4.  WAVE  APPROACH 

The  following  are  the  equations  that  govern  the  waves  in  the  layer  structure: 

V  x  V  x  Ej  -  kjEj  =  VjEj  j  =  !,■■■, N,  (6) 
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where  Vj  =  uj2fi<ij( r)  represents  the  volumetric  fluctuation  in  Region  j.  For  the  homogeneous  regions  above  and 
below  we  have 

V  x  V  x  E0  -  /cqEo  =  0,  (7a) 

V  x  V  x  Eat+i  —  k2N+l Ejv+i  =  0.  (7b) 

The  boundary  conditions  at  the  j-th  interface  are 

h  x  Ej(rj_,  Q)  =  hx  Ej+i(rj.,  Cj),  (8a) 

n  x  V  x  Ej  (rj_  ,  Q)  —  n  x  V  x  Ej+i(r_i_,  Q).  (8b) 

where  n  is  the  unit  vector  normal  to  the  j-th  interface  with  normal  pointing  into  the  medium  j.  This  system  is 
complemented  by  the  radiation  conditions  well  away  from  the  stack.  We  assume  that  we  know  the  solution  to 
the  problem  without  volumetric  fluctuations,  and  denote  it  as  E.  Let  the  Green’s  functions  to  this  problem  be 
denoted  as  Gy.  These  Green’s  functions  are  governed  by  the  following  set  of  equations: 

V  x  V  x  Gjfc(r,  r')  -  Gjk(r,  r')  =  lSjkS(r  -  r'),  (9a) 

n  x  Gjk(r±,Cj;r')  =  nx  G(j+1)k(r±,Q;r'),  (96) 

n  x  V  x  Gjk(rj_,Q;r')  =  h  x  V  x  G{j+1)k(r±,Q;r').  (9c) 

Another  pair  of  equations  similar  to  (9b)  and  (9c)  corresponding  to  (j  —  l)-th  interface  need  to  be  added  this 

list.  I  here  is  unit  dyad.  Using  these  Green’s  functions  and  the  radiation  conditions  the  wave  functions  can  be 
represented  as 


N  r 

Ei(r)  =  Ej(r)  +J2  dr'Gjk(r,  r>fc(r')Efc(r') 
fc= l 


(10) 


where  flk  =  {r^;  (k  <  z’  <  Cfc-i}-  Note  that  Vo  =  Vn+ i  =  0.  We  first  average  (10)  w.r.t.  volumetric  fluctuations 
to  get 


N  N 


(Ej(r))„  =  Ej(  +  [  dr1  I  dr"Gjk(r,  r')(GH(r',  r"))^(ufc(r,)ui(r"))(E/(r"))„.  (11) 

fc=i  i-i 

The  subscript  v  is  used  to  denote  averaging  with  respect  to  volumetric  fluctuations.  Here  we  have  used  a  first 
order  approximation  to  the  mass  operator  based  on  weak  fluctuations.  The  volumetric  fluctuations  in  different 
regions  are  assumed  to  be  uncorrelated,  which  means  that 

(vk(r')vl(r"))=5klCk(r'-r"),  (12) 

where  Ck  is  the  correlation  function  of  the  volumetric  fluctuations  in  Region  k.  Inserting  (12)  in  (11)  and 
employing  V  x  V  x  I  —  fe|  on  (11)  we  get 

V  x  V  x  (E,(r))„  -  kj(Ej(r))v  =  f  dr'(G^(r, t'))vCj(t  -  r')(EJ(r'))„.  (13) 

JQ,j 

Next  we  average  (13)  over  the  surface  fluctuations, 

V  x  V  x  (E,(r))^  -  kj{Ej(r))vs  =  f  dr’({ G,y  (r,  r'))„C,(r  -  r')(EJ(r'))1,) 

where  the  subscript  s  denotes  averaging  over  surface  fluctuations  and  =  {r'yj  Zj  <  z'  <  Zj- 1}.  We  approximate 
((Gjj(r, r'))vCj(r  -  r,)(Ei(r '))v^  as  (Gjj (r, r'))vsCj(r  -  r')(E_,-(r'))ra  and  obtain 


V  x  V  x  (E,(r))„a  -  k^ (Ej(r))ua  =  J  dr'(Gw(r,  r'))„sQ(r  -  r/)<Ej(r,))„a. 


(14) 
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We  call  this  the  weak  surface  correlation  approximation.  We  will  later  find  that  this  is  one  additional  approxima¬ 
tion  necessary  to  arrive  at  the  RT  system.  First  note  from  (14)  that  (Vx  Vx  I  —kj)(Ej(r))vs  =0forj  =  0,  iV+l. 
This  means  that  the  coherent  propagation  constants  in  regions  above  and  below  the  layer  stack  are  unaffected 
by  the  fluctuations  of  the  problem.  However,  they  indeed  get  modified  within  the  stack  region.  On  writing  (14) 

as  (V  x  Vx  I  —  k?  —  C)(ipj)  =  0,  where  C  denotes  the  integral  operator  f ^  dr'  (Gjj(r,  r'))vg  Cj( r  —  r'),  we  infer 
that  Xj  =  V(kj  +  £)  represents  the  mean  propagation  constant  in  j.  Observe  that  \j  depends  explicitly  on  the 
volumetric  fluctuations  in  Region  j  and  implicitly  on  the  fluctuations  of  the  stack,  both  volumetric  and  surface. 
This  is  in  contrast  to  the  RT  approach  where  7 j  depends  exclusively  on  the  volumetric  fluctuations  in  Region 
j.  Moreover,  \j  depends  on  the  polarization  if  the  fluctuations  of  the  problem  are  anisotropic.  Further,  even  if 
the  volumetric  fluctuations  are  isotropic  Xj  will  be  polarization  dependent  because  of  surface  reflections.  This 

is  in  contrast  to  the  RT  approach  where  1  j  is  polarization  dependent  only  when  the  volumetric  fluctuations  are 
anisotropic. 

Since  the  problem  is  statistically  homogeneous  in  azimuth  the  mean  fields  in  our  system  have  the  following 
form: 

(E^(r))«s  =  exp(zk_u  •  r)  jHP(k_u)p+  exp [iq?z\  +  Bp(k±i)pJ  exp[-zgpz]j  j  =  1, 2,  •  •  • ,  N,  (15) 

(Eg(r))„s  =  exp(iku  •  r)  jp^  exp [—ikoziz]  +  i?p(kj_j)po  exp[ik0ziz  ]},  (16) 

and 

(EAr+1(r))t>s  =  exp(ikj_j  •  r)Tp(k_Li)p“+1  exp [~ik(N+1)ziz\,  (17) 

where  the  superscript  p  stands  for  the  polarization,  either  horizontal  or  vertical,  p  is  the  unit  vector  representing 
polarization,  qj  is  the  z-coiriponent  of  Xj-  The  subscript  i  is  used  to  indicate  that  the  wave  vector  is  in  the 
incident  direction.  R  and  T  denote  respectively  the  mean  reflection  and  transmission  coefficient  of  the  stack.  Aj 
and  Bj  denote  respectively  the  mean  coefficients  of  up-going  and  down-going  waves  in  the  j-th  layer.  Based  on 
this  we  can  formulate  the  waves  averaged  w.r.t.  volumetric  fluctuations  as 


(Ej(r)>?  =  [  dk_Lexp(fkJ,-r)|H^(kj,  ,kj_i)q+exp[%z]+Hj9(k_L,k_Li)qi  exp[-i^z]|  j  =  1, 2,  •  •  • ,  N, 

(18) 

(E0(r))p  =  exp(ik_Lj  •  r)  exp[-ik0ziz}pQ  +  J  dk±  exp(?:k±  •  r)Rpq (k±,k±i)q^  exp [ik0zz\,  (19) 


and 

(E„+1(r»;  =  jij 

j  <ikj_  exp(?:k_L  •  r)TP9(kj_,  k_L?;)qAr+1  exp [~ik(N+i)zz\. 

(20) 

where  Aj, 
conditions 

Bj,  R  and  T  are  now  integral  operators  representing  scattering  from  rough  interfaces, 
associated  with  the  above  equations  at  the  j- th  interface  are 

The  boundary 

h  x 

(Ei  (r-L  tCj))v  —  x  (Ej+1  (r-L >Cj))v  j  —  1)  2,  •  •  jN 

(21a) 

and 

X 

<1 

X 

(Ej(r_L,  Cj))v  =fi  x  V  x  (Ey- 1  !  (r  Cy ) ) j  =  1,  2,  •  •  • ,  N. 

(216) 

The  above  system  may  be  solved  either  numerically  or  by  anyone  of  the  asymptotic  methods  available  in  rough 
surface  scattering  theory10-12  to  evaluate  the  mean  coefficients  that  appear  in  (15)-(17). 

We  proceed  now  to  the  analysis  of  the  second  moments,  by  starting  with  (10).  For  convenience  we  write  it 
in  symbolic  form  as 

N 

E,=E,+£g  jATfcEfc.  (22) 

fc= l 
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We  take  the  tensor  product  of  this  equation  with  its  complex  conjugate  and  average  w.r.t.  volumetric  fluctuations 
and  obtain 


N  N  N  N 

(Ej  0  E*)^  =  (Ej),  0  (E*)„  +  E  E  E  E<Gi*>"  ®  (Gjk')v^kk'W  (E,  ®  E; ?,)v  ,  (23) 

fc= 1  fc'= i  z=i  c=i 

where  K  is  the  intensity  operator  of  the  volumetric  fluctuations.  Employing  the  weak  fluctuation  approximation 
we  approximate  K  by  its  leading  term 


K kk>w  —  ( Vk  0  vl)Skk'w  I  •  (24) 

Next  we  average  (23)  w.r.t.  the  surface  fluctuations  and  employ  the  weak  surface  correlation  approximation  as 
before  to  get 


N 

(Ej  0  E*E  =  ((Ej),  0  (E*)„)  +  E  (  <Gjfc>„  0  (G*k)v  )  (vk  0  v*k)  (Efc  0  E%)va  .  (25) 

s  fc= l 

The  above  is  an  equation  for  the  second  moment  of  the  wave  function  E,  which  can  be  decomposed  into  a 
coherent  part  E  and  a  diffuse  part  E.  Therefore, 

(E  0  E*)  =  E  0  E*  +  (E  0  E*).  (26) 

The  coherent  part  is  not  of  much  interest;  we  know  that  it  is  specular  for  our  problem.  The  diffuse  or  the 
incoherent  part  is  of  more  interest.  Therefore  we  write  (25)  in  terms  of  diffuse  fields: 

N 

(Ej  0E‘)  =  ((Ej)w  0  (E^)  +  E  (  <Gjfc)„  0  <G *k)v  )  (vk  0  v*k)  (Efc  0  E*k)vs  ,  (27) 

S  1  S 

where  (E j)v  is  the  fluctuating  part  of  (E j)v.  Let  us  now  write  (27)  in  more  detail  as: 

(■Ej (r)  0  E*(r'))  =  ((E^r))„  0  (E*(?0)„)s+ 


+  E  /  dr,  I  dr'1((Gjfc(r,r1))tJ0(G*fe(r',r/1)>  \  \kk\4Ck(r,-r',))  (Efc(n)  0  E^ri)),,,  .  (28) 

k~ l  d Qfc 

As  it  stands  this  equation  is  very  difficult  to  solve  either  analytically  or  numerically.  Besides,  one  important  goal 
for  us  is  to  investigate  the  conditions  needed  for  employing  the  radiative  transfer  approach.  With  this  in  mind  we 
introduce  Wigner  transforms.  Note  that  (28)  is  an  equation  for  the  coherence  function  which  is  a  ‘space-space’ 
quantity.  On  the  other  hand  the  RT  equation,  as  we  saw  earlier,  is  an  equation  for  the  specific  intensity  which 
is  a  ‘phase-space’  quantity.  Wigner  transforms  serve  as  a  bridge  to  link  these  two  quantities. 18-21 

We  introduce  Wigner  transforms  of  waves  and  Green’s  functions  as 

£m  k )  =  y  d( r  -  r')  <Em(r)  0  E^r')}  (29) 


r  +  r 


2  ’ 


=  J  d(r-  r')  J  d(ri  —  r,1)e_ik^r_r ^e*1'*-ri_ri^(  (Gmn(r,  iq))^  0  (GJnn(r',  r^))^ 

(30) 


In  terms  of  these  transforms  (28)  becomes 


N 

£rn{ r,  k)  =  £^( r,  k)  +  E  \k™\4  J  dr'  J  da  J  d(3  C/mn(r,  k|r',  a)$n(a  -  /3)£„(r',  (3),  (31) 
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where  4>„  is  the  spectral  density  of  the  volumetric  fluctuations  in  the  n-th  layer.  We  have  used  the  superscript 
s  in  the  first  term  to  indicate  that  this  is  due  to  surface  scattering  as  defined  by  the  first  term  in  (28). 

The  fact  that  our  problem  has  translational  invariance  in  azimuth  implies  the  following: 


£m(r,k)  =  £m(z,  k), 

Gmn( r,  k|r',  1)  =  Gmn(z,  k| z' ,  l,  r_L  -  Tj_). 
Using  these  relations  in  (31)  we  have 

N 


£m{z,  k)  =  £sm(z,  k)  +  |fcn|4  J  dz1  J  da  J  d/3  Gmn(z,k\z',  a;  0)®n(a  - /3)£n(z',/3), 


(32a) 

(326) 


(33) 


where  Gmn{z,  k\z' ,  a;  0)  is  the  Fourier  transform  of  Gmn(z,  k^',  a;  rj_  —  r'±)  w.r.t.  rj_  —  evaluated  at  the 
origin  of  the  spectral  space.  To  proceed  further  we  need  to  evaluate  Gmn-  Further  we  need  to  relate  this  system 
with  that  of  RT,  which  involves  the  boundary  conditions  at  the  interfaces.  In  view  of  this  we  need  to  identify 
the  coherence  functions  corresponding  to  up-  and  down-going  wave  functions.  To  facilitate  this  we  decompose 
(■ Gmn)v  into  its  components, 


/  r-i  \  c  f^o  f^uu  /'-'i  ud  ,  s~^du  ,  dd 

\^mn/v  ~  °mn^rn  '  ^ mn  '  ^ ran  '  ^ mn  '  ^ mm 


(34) 


where  the  first  term  is  the  singular  part  of  the  Green’s  function.  The  superscripts  u  and  d  indicate  up-  and 
down-going  elements  of  the  waves.  The  other  components  are  due  to  reflections  from  boundaries.  These  are 
formally  constructed  using  the  concept  of  surface  scattering  operators  as  follows,12 


<G(t(  r,r')>r=  (2^)4  /  d j  dk'±{S%n(k±, 


(35) 


where  S“(n  is  the  surface  scattering  operator.  The  superscripts  a  and  b  on  S  are  used  to  indicate  whether  the 
waves  are  up-going  or  clown-going.  In  the  exponents,  a,  b  =  1  if  the  waves  are  up-going.  We  let  a,  b  =  —1  if 
the  waves  are  down-going.  The  ^-component  of  the  mean  propagation  constants  in  the  n-th  layer  is  denoted 
as  qn.  We  recall  that  Gmn  is  the  Wigner  transform  of  (  (Gmn)v  0  (G*nn)v  ^  .  The  superscripts  /i,  v  stand  for 

polarization,  either  h  or  v.  When  we  use  (34)  to  perform  the  Wigner  transform  we  ignore  all  cross  terms.  In 
other  words,  we  make  the  following  approximation, 


G 


mn 


—  SmnGm  +  G 


uu 

mn 


+  G 


ud 

mn 


+  G 


du 

mn 


+  G 


dd 

mni 


where  Q^n  is  the  Wigner  transform  of  ^  (G ^n)v  0  (G^*)^  y  .  This  approximation  reminds  us  of  the  concept 
of  incoherent  addition  of  intensities  associated  with  the  radiative  transfer  theory. 

With  the  introduction  of  this  representation  for  Gmn  hr  (31)  we  can  trace  up-  and  down-going  waves  to  obtain 
the  following  equations  for  the  coherence  function: 

C(2,k)  =  £™{z,k)  +  Jr^\km\A  J*  dz'  J  da  J  d/3  G>  (z,  k\z',  a,  0)$m(a  -  /3)£m(z',  /3) 


£dm{Zlk)  =  £^{z,k) 


N 


ifc”i4 


71=1 


da  J  d/3  GZn{zMz\<*^)$n{a-(3)£“{z',(3) 
d/3  G<(z,  k\z' ,  a ;  0)4>m(a  -  (3)£m{z' ,  (3) 
da  J  d/3G%n{z,k\z',a-,Vi)$>n{a-f3)£“{z',l3) 


(36  a) 


(366) 
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Note  that  summation  over  a  =  {w,  d}  is  implied  in  the  above  equations.  The  first  term  in  these  equations,  £sa, 
represents  the  contribution  due  exclusively  to  surface  scattering,  and  has  the  following  form: 


{C°(*,k)} 


His 


=  2t tS  \kz-  -a  [<  ( kj.)  +  q™( kjJ]  \  e 


Dia[^-  Qm]z 


c  ~  'i  mm  r  ~ 

{s^}  {s“:}  (kx,k±i) 


E  ,  ■  E*.  ■ 


(37) 

where  is  the  amplitude  of  the  up-going  wave  in  the  m-th  layer  after  volumetric  averaging  is  performed.  This 
means  that  it  is  a  random  function  of  surface  fluctuations.  When  we  substitute  (37)  and  the  expressions  for  Qmn 
in  (36)  we  find  that 


{£“(*, k)}^  =  2^8 \kz  -  Ia[C(kx)  +  C(kx)]}e“[<-^  {£“(*, kx)}^.  (38) 

On  substituting  this  in  (36)  and  differentiating  w.r.t.  z  we  obtain  the  following  transport  equations: 

+  *  fou(kx)  -  G£(kx)]  j  kJ-)  =^(^,kx)  + 


+  1^2  J  da-Ls»  |kX  -  ax;^  [^(kx)  +  ^(kx)]  -  ^a[<v(a-L)  +  Qt'(a±)]^£^vi{z,a±),  (39  a) 

{-^7  +  i  Mkx)  -  9C(k±)]|  k-L)  =  kx)+ 


\kn 


(2tt)5 


da±S* 


i 


kx  -  Mkx)  +5*(kx)]  -  ^a[q^(a±)  +  q*,{a±)}  [>  £^v,(z,  at±),  (39  b) 


where  £ represents  scattering  due  to  the  coherent  part  of  £,  whereas  the  integral  term  in  (39)  represents 

~a 

scattering  due  to  the  diffuse  part  of  £.  We  may  also  regard  £  as  the  source  to  our  transport  equations;  it  is 
given  as 


£fj,v  —  \kn 


T, 


kx  -  kXi;  [gM(kx)  +  g*(kx)]  -  ^ b  [gM(kx»)  +  <?*(kx*)] 


(S“®S“*)  :  [((S^)-E2))o((S^)-E ,»‘ 


(40) 


where  summation  over  b  is  implied.  Note  that  £a  in  (39)  includes  both  £u  and  £d  (corresponding  to  up  and 
down-going  waves).  When  the  superscripts  {a,  fo}  corresponds  to  u  the  value  of  {a,  b}  in  the  argument  of  4>m  take 
the  value  +1;  on  the  other  hand  when  the  superscripts  {a,  6}  corresponds  to  d  the  value  of  a  in  the  argument 
of  4>m  take  the  value  —1.  Since  all  quantities  in  (39)  and  (40)  correspond  to  the  same  layer  m  we  have  dropped 
the  subscript  m  to  avoid  cumbersome  notations.  To  obtain  appropriate  boundary  conditions  we  have  to  go  back 
to  the  integral  equation  representations  for  and  £dv  and  observe  their  behaviour  at  the  interfaces  and  try  to 
find  a  relation  between  them.  After  considerable  manipulations  we  managed  to  arrive  at  the  following  boundary 
conditions.  At  the  (m  —  l)-th  interface  we  have 


£i(Zm-n  kx)=  /  dkWftm-l.rr,  {^uKY)  ^m-uK) 


(41a) 


with  TZ=i{.  ®  R  where  Rm-i,m  is  the  stack  reflection  matrix  (not  the  local  reflection  matrix)  for  a  wave  incident 
from  below  on  the  (m  —  l)-th  interface.  Similarly 


kx) 


(kx.k'x))^^.! 


k'x), 


(416) 


where  IZm+i  ,m  is  the  tensor  product  of  stack  reflection  matrix  for  a  wave  incident  from  above  on  the  m-th 
interface.  We  were  able  to  obtain  the  boundary  conditions  only  after  imposing  certain  approximations  such  as 
the  one  given  below.  Consider  the  following  identity: 

s^m  =  Fm  Rrn-i.m  (S>  +  srm}  Fm  (42) 
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where  Fm  =  diag  {elqhdrn ,elqvdm}.  Notice  that  this  is  an  operator  relation  where  all  elements  are  operators. 
Taking  the  tensor  product  of  (42)  with  its  complex  conjugate  we  have 


G*du  cdu* 
^mm  ^  ^mm 


(Fr 


i)  (Rm— l,r; 


R 


'771—1 


,m)  ({^ 


l}®{S>+Srm}1(Fm®F^). 


Next  we  average  (43)  w.r.t.  surface  fluctuations  and  get 


(43) 


(s%m  0  *  (Fm  0  F*J  (Rm_1>ro  0  Rl-i,m)  (  {S>  +  SZJ  ®  {S>  +  S““ro}*  ^  (Fm  0  F*m)  (44 a) 


where  we  have  approximated  that  the  two  tensor  products  in  the  middle  are  weakly  correlated.  A  further 
approximation  that  we  impose  is  given  as  follows 


{s>  +  srm}  ®  {s>  +  srm}  )  *  s>  0  s>*  +  ( srm  ®  sn 


These  are  the  kinds  of  approximations  required  to  arrive  at  our  boundary  conditions. 


(446) 


5.  TRANSITION  TO  RADIATIVE  TRANSFER 

Now  we  have  to  transition  from  this  transport  equation  (44)  to  the  phenomenological  radiative  transfer  equation 
discussed  earlier.  To  accomplish  this  we  have  to  link  the  key  quantities  of  waves  and  radiative  transfer,  viz., 
coherence  function  and  specific  intensity.  The  relation  between  them  is  obtained  by  computing  the  energy  density 
using  the  two  concepts.  Thus  we  have 


ie(EM(r)E*(r))  =  JjTt  J  cttV^r,  s). 

Wigner  transform  provides  us  following  relation: 

(E^r)Et(r))  =  ^  J  dk±£^(z,k±). 
From  (45)  and  (46)  we  relate  I  to  £  as 


s) 


1  k'2 
2rj  (2-7t)2 


cos  9£liV(z,k±). 


(45) 


(46) 


(47) 


There  is  still  one  difference  that  needs  to  be  ironed  out  before  we  transition  to  the  RT  equations.  Notice  that  in 
our  wave  approach  we  obtained  transport  equations  for  6,  which  is  the  fluctuating  part  of  the  coherence  function. 
On  the  other  hand  the  phenomenological  RT  equations  are  traditionally  written  for  total  intensities.  Therefore 
we  have  to  express  our  transport  equations  in  terms  of  £.  Notice  that  £  =  £  +  £,  where  £ ,  the  average  part  of 
£  satisfies: 

{ —  ia(q^  ~  it) }  kj_)  =  °i  (48) 

Using  (48)  in  (39)  we  obtain 

j^  +  iMk-L)  -^(kj-)]}^(^kJ-)  =  |^|2  J  da±S>®S>* 

jk_L  -  ar,  ^  [9M(kj.)  -  g*(kj_)]  -  ^a[qli'(a±)  -  g*/(aj_)]|  £^v>(z,  a±),  (49 a) 

+i[gAt(kj-)  -^(kJ-)]}^(2,kj_)  =  |^2  J  da±S<®S<* 

4>m  |k±  -  aii  ~  Mk±)  -  qt( k±)]  -  ^fl[v(“i)  -  j  ££>„>&,  c*_l),  (496) 
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Notice  that  this  equation  is  expressed  entirely  in  total  intensity.  Now  we  can  transition  to  the  phenomenological 
RT  equations.  Using  the  relation  between  £  and  I  we  change  the  integration  variable  to  solid  angle  and  arrive 
at  the  following  equation, 


jcos^  +  rfej  If(z,s)  =  j  dn'P^(Q,Q')Ij(z,s'), 
j— COS^+7 lij}lf(Z,s)  =  I 


(50a) 

(506) 


where  V  is  the  extinction  matrix  and  P  is  the  phase  matrix.  Implicit  summation  over  superscripts  a  and  subscript 
v  is  assumed  in  (50).  To  facilitate  comparison  with  the  results  of  Ulaby  et  al., 6  and  Lam  and  Ishimaru'  we  have 
used  a  modified  version  of  Stokes  vector.17  Instead  of  the  standard  form  {/,  Q,  U,V}  we  use  {(/  +  Q)/ 2,  (/  — 
<3)/2,  U,  V}.  The  subscript  of  I  denote  the  element  number  of  our  Stokes  vector.  Although  the  structure  of  this 
equation  is  identical  to  that  of  RT  (equation  (2))  the  elements  of  the  phase  matrix  and  the  extinction  matrices 
are  not  the  same.  The  primary  reason  is  because  of  the  differences  in  the  real  part  of  the  mean  propagation 
constants  of  horizontally  and  vertically  polarized  waves.  On  assuming  that  q'h  =  q'v  =  k'mz  we  obtain  the  following 
expressions  for  the  extinction  and  phase  matrices: 


??=  -  cos  6  diag  {2 q”,  2q'^  q "  +  q%,  q "  +  q£}  (51a) 

p$  =  \\km\^m  {k±  -  k'±;  k'm [a cos 9  6cos*']}^  (516) 

For  n  =  v,  h, 

Kt  =  (pa-v'b)2  Kh  =  {pa-  h'6)2 

Ku  =  (/*“  ■  K  (m°  •  h'b)  V%  =  0  (52) 

Similarly, 

Kv  =  2  (va  •  Vb)  (ha  •  Vb )  Vtfh  =  2  (va  •  h'6)  (ha  •  h'b) 

Vtfu  =  (va  •  Vb)  (ha  •  h'&)  +  (va  •  h'6)  (h°  •  v'&)  V?jbv  =  0  (53) 

nbv  =  Kh  =  Ku  =  o 

v^v  =  (va  •  v,b)  (h°  •  h'&)  -  (va  •  h,b)  (ha  •  v,&)  (54) 

Noting  the  implied  summation  over  a  in  (50)  we  see  that  they  are  identical  to  the  RT  equations  given  in  Section  2. 
Now  we  have  explicit  expressions  for  the  extinction  matrix  and  phase  matrix  in  terms  of  the  statistical  parameters 
of  the  problem  thanks  to  our  wave  approach.  We  next  turn  our  attention  to  the  boundary  conditions.  In  our 
wave  approach  we  obtained  BCs  in  terms  of  ‘stack’  reflection  matrix  R,  whereas  in  the  RT  approach  the  BCs  are 
given  in  terms  of  the  local  interface  reflection  matrices.  We  can  readily  reconcile  with  this  apparent  difference. 
Note  that  the  BC  in  the  wave  approach  forms  a  closed  system  whereas  in  the  RT  approach  it  is  open  (linked  to 
adjacent  layer  intensities).  Let  us  take  a  look  at  the  BC  at  the  (m  —  l)-th  interface.  Rm-2,m  can  be  expressed 
in  terms  of  Rm_2,m,-i  as  follows, 


Rm— l,m —  R 


■m—  l,m 


-1 


-l,m* 


(55) 


This  is  the  relation  between  the  stack  reflection  coefficients  of  adjacent  interfaces.  The  R  and  T  are  local  (single 
interface)  reflection  and  transmission  matrices  at  the  (m  —  l)-th  interface.  On  operating  Ejj,  with  (54)  we  get 

=  kra-l,mE“  +  Tm>m_iE  ^J_1.  (56) 

Notice  that  this  boundary  condition  now  involves  only  local  interface  Fresnel  coefficients. 
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(57) 


Similarly  we  write  Rm+i,-m  in  terms  of  Rm+2,m+i  and  hence  obtain  the  BC  at  the  m-th  interface 

EU  D  T 7id  |  rri  TTU 

m  m,m+l -^m+l  * 

Take  the  tensor  product  of  (56)  with  its  complex  conjugate  and  average  w.r.t.  surface  fluctuations.  Employing 
the  Wigner  transform  operator,  we  obtain  boundary  condition  at  the  (m  —  l)-th  interface  similar  to  that  of  the 
RT  system.  However,  the  reflection  and  transmission  matrices  used  in  the  RT  system  correspond  to  unperturbed 
medium  as  opposed  to  the  average  medium  as  in  the  case  of  the  wave  approach. 

6.  DISCUSSION 

Now  that  we  have  made  the  transition  from  statistical  wave  theory  to  radiative  transfer  theory  it  is  instructive 
to  itemize  the  assumptions  implicitly  involved  in  the  RT  approach. 

1.  Weak  fluctuations 

2.  Quasi-stationary  field  approximation. 

3.  Statistical  homogeneity  of  fluctuations. 

These  are  the  three  well-known  conditions  necessary  for  the  unbounded  random  medium  problem.  However,  if 
the  medium  is  bounded  we  need  to  impose  additional  conditions.  We  found  that  the  extinction  coefficients  cal¬ 
culated  in  the  wave  approach  and  the  RT  approach  are  different  and  only  after  applying  further  approximations 
can  they  be  made  to  agree  with  each  other.  The  following  two  additional  conditions  are  required  for  bounded 
random  media: 

4.  While  calculating  effective  propagation  constants  we  ignore  terms  involving  interface  reflection  coefficients. 

In  other  words  we  ignore  all  boundary  effects. 

5.  While  carrying  out  Wigner  transforms  of  the  Green’s  functions  we  need  to  neglect  the  interactions  between 
components. 

6.  Use  unperturbed  medium  parameters  while  imposing  boundary  conditions. 

The  above  two  conditions  are  required  when  the  boundaries  are  planar.  However,  if  they  are  randomly  rough  we 
need  to  further  impose  two  more  conditions: 

7.  Weak  surface  correlation  approximation. 

8.  All  fluctuations  of  the  problem  are  statistically  independent. 

To  summarize,  we  have  enquired  into  the  assumptions  involved  in  adopting  the  radiative  transfer  approach 
to  scattering  from  layered  random  media  with  rough  interfaces.  To  facilitate  this  enquiry  we  adopted  a  wave 
approach  to  this  problem  and  derived  the  governing  equations  for  the  first  and  second  moments  of  the  wave 
fields.  We  employed  Wigner  transforms  and  transitioned  to  the  system  corresponding  to  that  of  radiative  transfer 
approach.  In  this  process  we  found  that  there  are  more  conditions  implicitly  involved  in  the  RT  approach  to  this 
problem  than  it  widely  believed  to  be  necessary.  With  the  recent  development  of  fast  and  efficient  algorithms 
for  scattering  computations  and  the  enormous  increase  in  computer  resources  it  now  feasible  to  take  an  entirely 
numerical  approach  to  this  problem  without  imposing  any  approximations.  However,  to  keep  the  size  of  the 
problem  to  be  manageable  only  special  cases  have  been  studied  thus  far. 22-24  Hence  it  is  very  much  of  relevance 
and  interest  to  use  the  RT  approach  to  these  problems.  However,  one  should  keep  in  mind  the  assumptions 
involved  in  such  an  approach. 
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